knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(message = FALSE)
knitr::opts_chunk$set(warning = FALSE)
knitr::opts_chunk$set(cache = TRUE)
library(tidyverse)
library(cowplot)
library(broom)
library(readxl)
library(janitor)
library(plotrix)
library(here)
library(growthTools)
library(rootSolve)
library(nlstools)
library(nls.multstart)
source(here("R-scripts", "predict_monod.R"))
Read in data
treatments <- read_excel(here("data-general", "ChlamEE_Treatments_JB.xlsx")) %>%
clean_names() %>%
mutate(treatment = ifelse(is.na(treatment), "none", treatment)) %>%
filter(population != "cc1629")
nitrate <- read_csv(here("data-processed", "nitrate-abundances-processed.csv")) %>%
group_by(population, nitrate_concentration, well_plate) %>%
mutate(N0 = RFU[[1]]) %>%
filter(population != "COMBO")
# Eyeball exponential 2 day cutoff ----------------------------------------
nitrate_exp <- nitrate %>%
mutate(exponential = case_when(days < 2 ~ "yes",
TRUE ~ "no")) %>%
filter(exponential == "yes")
# nitrate_exp %>%
# ggplot(aes(x = days, y = RFU, color = factor(nitrate_concentration))) + geom_point() +
# facet_wrap( nitrate_concentration~ population, scales = "free_y") + scale_color_viridis_d()
# ggsave(here("figures", "nitrate-exp.png"), width = 45, height = 45)
growth_rates <- nitrate_exp %>%
group_by(nitrate_concentration, population, well_plate) %>%
do(tidy(nls(RFU ~ N0 * exp(r*days),
data= ., start=list(r=0.01),
control = nls.control(maxiter=100, minFactor=1/204800000)))) %>%
ungroup()
Fit the Monod model to the growth rates
monod_fits <- growth_rates %>%
mutate(nitrate_concentration = as.numeric(nitrate_concentration)) %>%
group_by(population) %>%
do(tidy(nls(estimate ~ umax* (nitrate_concentration / (ks+ nitrate_concentration)),
data= ., start=list(ks = 1, umax = 1), algorithm="port", lower=list(c=0.01, d=0),
control = nls.control(maxiter=500, minFactor=1/204800000))))
bs_split <- monod_fits %>%
select(population, term, estimate) %>%
dplyr::ungroup() %>%
spread(key = term, value = estimate) %>%
split(.$population)
all_preds_n <- bs_split %>% ### here we just use the fitted parameters from the Monod to get the predicted values
map_df(predict_monod, .id = "population")
all_predsn_2 <- left_join(all_preds_n, treatments, by = c("population")) %>%
distinct(ancestor_id, treatment, nitrate_concentration.x, .keep_all = TRUE)
all_growth_n2 <- left_join(growth_rates, treatments) ## changed this to the 'cut' version. can switch back later
all_growth_n2 %>%
# mutate(estimate = mu) %>%
mutate(nitrate_concentration = as.numeric(nitrate_concentration)) %>%
# filter(treatment == "N", ancestor_id == "anc3") %>%
ggplot(aes(x= nitrate_concentration, y= estimate)) +
geom_point() +
# geom_point(aes(x = nitrate_concentration.x, y = estimate), data = filter(nitrate_eyeball_exp, population == 27), color = "blue") +
geom_line(data=all_predsn_2, aes(x=nitrate_concentration.x, y=growth_rate, color = treatment), size = 1) +
facet_grid(treatment ~ ancestor_id) +
geom_hline(yintercept = 0.1, linetype = "dotted") +
ylab("Exponential growth rate (/day)") + xlab("Nitrate concentration (uM)")
ggsave(here("figures", "nitrate-monod-curves-eyeball-2day.png"), width = 12, height =8)
monod_wide <- monod_fits %>%
select(population, term, estimate) %>%
spread(key = term, value = estimate)
m <- 0.1 ## set mortality rate, which we use in the rstar_solve
Find R*
rstars <- monod_wide %>%
mutate(rstar_solve = ks*m/(umax-m)) ## analytical
rstars2 <- left_join(rstars, treatments, by = "population") %>%
distinct(population, ks, umax, .keep_all = TRUE)
rstars2 %>%
group_by(treatment) %>%
summarise_each(funs(mean, std.error), rstar_solve) %>%
ggplot(aes(x = reorder(treatment, mean), y = mean)) + geom_point() +
geom_errorbar(aes(ymin = mean - std.error, ymax = mean + std.error),width = 0.1) +
ylab("R* (umol N)") + xlab("Selection treatment") + geom_point(aes(x = reorder(treatment, rstar_solve), y = rstar_solve, color = ancestor_id), size = 2, data = rstars2, alpha = 0.5) +
scale_color_discrete(name = "Ancestor")
ggsave(here("figures", "nitrate-r-star-exponential-2day.png"), width = 6, height = 4)
# sequential exponential fits ---------------------------------------------
fitting_window <- function(x) {
growth_rates <- nitrate %>%
top_n(n = -x, wt = days) %>%
group_by(nitrate_concentration, population, well_plate) %>%
do(tidy(nls(RFU ~ N0 * exp(r*days),
data= ., start=list(r=0.01),
control = nls.control(maxiter=100, minFactor=1/204800000)))) %>%
mutate(number_of_points = x) %>%
ungroup()
}
windows <- seq(4,7, by = 1)
multi_fits <- windows %>%
map_df(fitting_window, .id = "iteration")
multi_fits %>%
ggplot(aes(x = number_of_points, y = estimate, group = well_plate)) + geom_point() + geom_line() +
facet_wrap( ~ nitrate_concentration)
exp_fits_top <- multi_fits %>%
group_by(well_plate) %>%
top_n(n = 1, wt = estimate)
Fit the Monod model to the growth rates
monod_fits_seq <- exp_fits_top %>%
mutate(nitrate_concentration = as.numeric(nitrate_concentration)) %>%
group_by(population) %>%
do(tidy(nls(estimate ~ umax* (nitrate_concentration / (ks+ nitrate_concentration)),
data= ., start=list(ks = 1, umax = 1), algorithm="port", lower=list(c=0.01, d=0),
control = nls.control(maxiter=500, minFactor=1/204800000))))
bs_split_seq <- monod_fits_seq %>%
select(population, term, estimate) %>%
dplyr::ungroup() %>%
spread(key = term, value = estimate) %>%
split(.$population)
all_preds_seq <- bs_split_seq %>% ### here we just use the fitted parameters from the Monod to get the predicted values
map_df(predict_monod, .id = "population")
all_preds_seq2 <- left_join(all_preds_seq, treatments, by = c("population")) %>%
distinct(ancestor_id, treatment, nitrate_concentration.x, .keep_all = TRUE)
all_growth_seq2 <- left_join(exp_fits_top, treatments) ## changed this to the 'cut' version. can switch back later
all_growth_seq2 %>%
mutate(nitrate_concentration = as.numeric(nitrate_concentration)) %>%
ggplot(aes(x= nitrate_concentration, y= estimate)) +
geom_point() +
geom_line(data=all_preds_seq2, aes(x=nitrate_concentration.x, y=growth_rate, color = treatment), size = 1) +
facet_grid(treatment ~ ancestor_id) +
geom_hline(yintercept = 0.1, linetype = "dotted") +
ylab("Exponential growth rate (/day)") + xlab("Nitrate concentration (uM)")
ggsave(here("figures", "nitrate-monod-curves-sequential-method.png"), width = 12, height =8)
monod_wide_seq <- monod_fits_seq %>%
select(population, term, estimate) %>%
spread(key = term, value = estimate)
m <- 0.1 ## set mortality rate, which we use in the rstar_solve
Find R*
rstars_seq <- monod_wide_seq %>%
mutate(rstar_solve = ks*m/(umax-m)) ## analytical
rstars2_seq <- left_join(rstars_seq, treatments, by = "population") %>%
distinct(population, ks, umax, .keep_all = TRUE)
rstars2_seq %>%
group_by(treatment) %>%
summarise_each(funs(mean, std.error), rstar_solve) %>%
ggplot(aes(x = reorder(treatment, mean), y = mean)) + geom_point() +
geom_errorbar(aes(ymin = mean - std.error, ymax = mean + std.error),width = 0.1) +
ylab("R* (umol N)") + xlab("Selection treatment") + geom_point(aes(x = reorder(treatment, rstar_solve), y = rstar_solve, color = ancestor_id), size = 2, data = rstars2_seq, alpha = 0.5) +
scale_color_discrete(name = "Ancestor")
ggsave(here("figures", "nitrate-r-star-sequential-method.png"), width = 6, height = 4)
# logistic ----------------------------------------------------------------
fits_many_n0 <- nitrate %>%
group_by(population, well_plate) %>%
nest() %>%
mutate(fit = purrr::map(data, ~ nls_multstart(RFU ~ K/(1 + (K/N0- 1)*exp(-r*days)),
data = .x,
iter = 500,
start_lower = c(K = 100, r = 0),
start_upper = c(K = 500, r = 2),
supp_errors = 'N',
na.action = na.omit,
lower = c(K = 10, r = 0),
upper = c(K = 4000, r = 4),
control = nls.control(maxiter=1000, minFactor=1/204800000))))
fits_many <- fits_many_n0
fits_well_plate <- fits_many %>%
select(well_plate)
key <- nitrate %>%
select(well_plate, population, nitrate_concentration) %>%
distinct(well_plate, .keep_all = TRUE)
pops <- left_join(fits_well_plate, key)
info <- fits_many %>%
unnest(fit %>% map(glance))
# get params
params <- fits_many %>%
filter(fit != "NULL") %>%
unnest(fit %>% map(tidy))
# write_csv(params, here("data-processed", "logistic-parameters-nitrate.csv"))
bad_fits <- params %>%
filter(term == "K", estimate == 4000)
bad_fits2 <- left_join(bad_fits, key)
# CI <- fits_many %>%
# filter(fit != "NULL") %>%
# unnest(fit %>% map(~ confint2(.x) %>%
# data.frame() %>%
# rename(., conf.low = X2.5.., conf.high = X97.5..))) %>%
# group_by(., well_plate) %>%
# mutate(., term = c('r', 'K')) %>%
# ungroup()
#
# # merge parameters and CI estimates
# params <- merge(params, CI, by = intersect(names(params), names(CI)))
new_preds <- nitrate %>%
do(., data.frame(days = seq(min(.$days), max(.$days), length.out = 15), stringsAsFactors = FALSE))
# max and min for each curve
max_min <- group_by(nitrate, well_plate) %>%
summarise(., min_days = min(days), max_days = max(days)) %>%
ungroup()
# create new predictions
# preds2b <- fits_many %>%
# unnest(fit %>% map(augment, newdata = new_preds)) %>%
# merge(., max_min, by = 'well_plate') %>%
# group_by(., well_plate) %>%
# filter(., days > unique(min_days) & days < unique(max_days)) %>%
# rename(., RFU = .fitted) %>%
# ungroup()
preds <- fits_many %>%
unnest(fit %>% map(augment))
new_preds <- nitrate %>%
do(., data.frame(days = seq(min(.$days), max(.$days), length.out = 15), stringsAsFactors = FALSE))
# max and min for each curve
max_min <- group_by(nitrate, well_plate) %>%
summarise(., min_days = min(days), max_days = max(days)) %>%
ungroup()
key <- nitrate %>%
select(well_plate, population, nitrate_concentration) %>%
distinct(well_plate, .keep_all = TRUE)
preds3 <- left_join(preds, key, by = c("well_plate", "population")) %>%
left_join(., treatments) %>%
distinct(.fitted, population, nitrate_concentration, well_plate, .keep_all = TRUE)
# preds3 <- left_join(preds2, key, by = c("well_plate", "population"))
nitrate2 <- left_join(nitrate, treatments)
ggplot() +
geom_point(aes(days, RFU, color = factor(nitrate_concentration)), size = 2, data = nitrate2) +
geom_line(aes(days, .fitted, group = well_plate, color = factor(nitrate_concentration)), data = preds3) +
facet_grid(ancestor_id ~ treatment, scales = "free") +
scale_color_viridis_d() +
ylab('RFU') +
xlab('Days')